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We study spin diffusion and spin waves in paramagnetic quantum crystals 
(solid 3 He, for example) by direct simulation of a square lattice of atoms in- 
teracting via a nearest-neighbor Heisenberg exchange Hamiltonian. Recently, 
Cowan and Mullin have used a moments method to study spin transport at 
arbitrary polarizations. We test their analytic results by calculating the sta- 
tistical spin correlation function from molecular dynamics simulations using 
a Monte Carlo algorithm to average over initial spin configurations. Since it 
is not practical to diagonalize the S = ^ exchange Hamiltonian for a lattice 
which is of sufficient size to study long-wavelength (hydrodynamic) fluctua- 
tions, we instead study the S — > oo limit and treat each spin as a vector with 
a classical equation of motion. We compare our simulations with the as- 
sumptions of the moments method regarding the short-time behavior and the 
long exponential tail of the correlation function. We also present our numer- 
ical results for the polarization dependence of the longitudinal spin diffusion 
coefficient and the complex transverse spin diffusion coefficient. 



1. INTRODUCTION 

The Leggett-Rice (LR) equation^ was originally derived for paramag- 
netic Fermi fluids, but it has been shown that it applies to the spin dynam- 
ics of any polarized system interacting via quantum exchange. 2 In the lin- 
ear regime, the LR equation describes conventional diffusion of longitudinal 
Fourier modes S z oc e ici ' r ~ D ^ q *, where Dn is the longitudinal diffusion coeffi- 
cient, as well as damped transverse spin waves S + = S x + iS y oc e *q r - £)+ '? 2i ) 
characterized by a complex diffusion coefficient D + . The transverse diffusion 
coefficient can be written in terms of a bare diffusion coefficient as 
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where D± is real and the denominator takes into account the precession 
of the spin current about the molecular field with quality factor fiM = 
lmD + /ReD + . 

Cowan and Mullin et aZpH21 (CM) have applied this description to the 
spin dynamics of paramagnetic 3 He in its hep phase, with a nearest-neighbor 
pair exchange Hamiltonian, and have used a moments method, originally due 
to RedfielcP and deGennesj^to calculate the polarization dependence of the 
transport coefficients. At zero polarization P = they found ordinary spin 
diffusion (fxM = 0, D± = Dm ~ uo 2 ), where ui is the exchange frequency and 
a is the lattice spacing. For P — > 1 they found that all the transport coeffi- 
cients diverge, but with D±//iM finite, as one would expect for undamped 
spin waves. Over a wide range of P they found D± ~ Dm , although the two 
diverge differently as P — ► 1. 

In order to test this picture we carried out a numerical study of spin 
diffusion on a 2d square lattice with the Hamiltonian 

H = -Y,BSl-Y,JSi-Sj (2) 

i i,j 

were the second sum is over nearest neighbors. We shall restrict our attention 
to the weakly paramagnetic case, with J/T — > but B/T remaining finite. 
Since a lattice of N spins involves 2^ quantum states, it is not practical 
to diagonalize the S = \ exchange Hamiltonian for a lattice which is of 
sufficient size to study long-wavelength (hydrodynamic) fluctuations. We 
instead study the S — > oo limit and treat each spin as a vector S(r) with a 
classical equation of motion, 

|s(r)=u,£S(r)xS(r') (3) 

r' 

where the sum is over nearest neighbors. For convenience, we have moved 
to the Larmor frame, and we set the exchange frequency oj = J/% and the 
vector magnitude of the spins |S| equal to 1. The longitudinal diffusion 
coefficient was calculated by simulating the motion of a hydrodynamic {q — > 
0) fluctuation Sg(t). We define the correlation function^ 

G z q (t) = (Sl(t)S z q *m (4) 

which is assumed to have an exponential decay oc e- D W q 1 as t -» oo. The 
brackets indicate an average over initial conditions sampled from an equilib- 
rium ensemble. 
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Fig. 1. Log plot of the longitudinal spin correlation function G z q {t) at zero 
polarization. The slope is equal to —D\\q 2 from which a value of Dm = 
0.51 ±0.01 (in units of o;a 2 )is obtained. To the right, the short time behaviors 
of the first and second derivatives are also shown. For comparison, the 
dashed line in the G" plot is the Gaussian ansatz used by CM. 3 



2. NUMERICAL SIMULATIONS 

We shall outline the procedure for calculating Dm. First, an initial 
spin configuration was obtained by generating a lattice of randomly oriented 
spins Sj with a Boltzman distribution e^ B ' T ^ s i . The exchange energy was 
neglected in the Boltzman factor, since in a weakly paramagnetic system 
the spins are basically uncorrelated, although their motion is not. The spin 
configuration was then allowed to evolve according to Eq|3J and the spin 
correlation function G z q {t) was calculated from the definition Eq0J The 
average was then performed over N initial configurations were N = 1000 
was typical. The size of the lattice was L x x L y = 64 x 8 spins with periodic 
boundary conditions. The length L x of the lattice was chosen so that the 
wave number q = 2n/L x = 0.098 was small. To insure that the lattice was 
a 2-d system, L y was chosen so that the mean free time r = Dy/fwa) 2 was 
shorter than the time L 2 /(D||7r 2 ) to diffuse across the width of the lattice. 
The equation of motion was integrated with a second-order modified Euler 
method, with time step At = 0.01. The typical motion of a spin in a time 
step was half a degree or so. The conserved quantities such as the total 
energy, the polarization, and the individual spin magnitudes were found to 
be constant to a high degree of accuracy. 

In Fig.l the correlation function G z q {t) is shown for P = 0. A similar 
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Fig. 2. Plot of log|£r+(i)| (solid line)and the phase (dashed line) from the 
previous figure. A quality factor of \iM = 4.6 ±0.1 is found from the ratio 
of the slopes of the two lines. The inset shows the instantaneous value of 
fiM. 

procedure was used to calculate D± and the quality factor [iM from 5+ 
Some typical results are shown in Fig. 2 for P = 0.673. 



3. RESULTS AND DISCUSSION 



The polarization dependence of the transport coefficients are shown in 
Fig. 3-4. The results are very similar to results of CM obtained with a mo- 
ments method. As expected, all the transport coefficients diverge as P — ► 1, 
although D±/fiM remains finite. There was no indication of anisotropy in 
the diffusion coefficients for P < 0.8. At higher polarizations the divergence 
of -Dy required prohibitively large lattices and calculation times. This was 
not a problem in the transverse case, since ReZ) + — > as P — > 1. Some 
disagreement with CM is expected since they considered a 3d spin- 1/2 sys- 
tem using a mean field approximation. When our transport coefficients are 
plotted vs. 1 — P 2 on a log-log plot the resulting curve fits are linear and 
indicate a divergence cx (1 — p 2 ) 3 0±0 - 3 . This disagrees with CM who found 
an exponent of —1 for the transverse coefficients and —1/2 for -Dy. 5 

However, there is remarkable agreement at low polarization. For ex- 
ample, at zero polarization CM obtained Dn = D± = 0.79 for a 2d square 
latticed This is higher than our value of 0.51, but when comparing classical 
and quantum results one should make the substitution S 2 — > S(S + 1) with 
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Fig. 3. Polarization dependence of D\\ (squares) and Z?j_ (triangles). The 
error bars are smaller than the size of the symbols, except for the two D» 
points at P > 0.8. The inset shows Du and the effective transverse diffusion 
coefficient KeD + for comparison with Fig.l of RefEJ 
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Fig. 4. Polarization dependence of the quality factor. The inset shows ImD + . 
Obviously, ImD + ~ Pujcl 2 . Since the spin wave frequency is lmD + q 2 , this 
supports the notion of a molecular precessional field. It also agrees with the 
mean field results of CM, and agrees with the exact solution at P = 1. 
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S = 1/2. We omit the details of the analysis here, but this increases our 
value by a factor of \/3 to 0.87. Also, for 3d hep they obtain fj,M/P = 2.3 
as P — > 0, 5 which is higher than our value of 1.5. The quantum correction 
above cancels for this quantity, but from connectivity arguments, one would 
expect the molecular field in a 3d system to be stronger by a factor 3/2, 
which would increase our value to 2.3, in perfect agreement with CM. 

Finally, our simulations with classical spin vectors support some of the 
assumptions of the moments method. For long times, hydrodynamic fluctu- 
ations were found to decay exponentially. For short times, the second deriva- 
tive of the correlation function was well-described by a Gaussian ansatz and 
decayed on a time scale r = D^/(uja) 2 (See Fig.l). 
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